home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / procssng / ccs / ccs-11tl.lha / lbl / hips / sources / mahe / ahecalc.c < prev    next >
Encoding:
C/C++ Source or Header  |  1991-09-17  |  10.6 KB  |  406 lines
  1. /* 5/30/1988 Modified by Shie-Jue Lee:                */
  2. /*    1. data type for equal_const                */
  3. /*    2. Iteration controls for histogram calculation        */
  4.  
  5. typedef long GREYTYPE;
  6.  
  7. #include <stdio.h>
  8. #include <signal.h>
  9. #include <math.h>
  10. #include "ahe.h"
  11.  
  12. #define Malloc(x) malloc((unsigned)(x))
  13.  
  14. #define iptr1(Y,X)    (im1+((Y) * numline)+(X))
  15. #define iptr2(Y,X)    (im2+((Y) * numline)+(X))
  16.  
  17. #define mapval(J,K)    (*(maps + ((J) * regline) + (K)))
  18.  
  19. int         status;        /* progress of calc. indicators     */
  20.  
  21. ahecalc(im1, im2, dimv, minmax, nreg, clipfrac, argv)
  22.      /* CHANGED */
  23.      unsigned char *im1,            /* input image     */
  24.      *im2;            /* output image */
  25.      int *dimv,            /* array of image bounds     */
  26.      minmax[2],        /* minimum and maximum inputs     */
  27.      nreg[DIMNO];        /* number of regions in x, y, z  */
  28.      float clipfrac;        /* fraction of average bucket cnt */
  29.      char **argv;        /* command line for status report via 'ps' */
  30. {
  31.  
  32.  
  33. /* CHANGED */
  34.     int tmpint ;
  35.     int         min,        /* minimum grey value             */
  36.                 max,        /* maximum grey value             */
  37.                 maxlevels,    /* number of grey levels in image     */
  38.                 numregions,    /* number of regions to be used         */
  39.                *del[DIMNO],    /* region sizes in each dimension     */
  40.                *curdel,        /* current region under consideration     */
  41.                 regsize,    /* compute size of current region     */
  42.                 excess,        /* excess pixels in each dimension     */
  43.                 newdel[DIMNO],    /* size of current regions         */
  44.                 start[DIMNO],    /* current starting position in image     */
  45.                 regline,    /* number of regions in line of image     */
  46.                 cliplimit,    /* computed limit of bucket count     */
  47.                 limit,        /* cliplimit - ptemp             */
  48.                 numline,    /* number of pixels per line of image     */
  49.                *histo,        /* histogram array             */
  50.                 ptemp;        /* floor of p                  */
  51.     float      *cumhisto,    /* cummulative histogram array           */
  52.                 p;        /* base for cumulative histogram     */
  53.  
  54.     /* CHANGED */
  55.     register unsigned char *i1,    /* fast pointers to subsegments of      */
  56.                *i2;        /* image arrays                 */
  57.  
  58.     int         ymin,        /* bounds for mapping of image         */
  59.                 ymax,
  60.                 xmin,
  61.                 xmax;
  62.  
  63.     int         jval,        /* map indicies                 */
  64.                 kval;
  65.     double      recy,        /* constants to step thru image in mapping */
  66.                 recx,
  67.                 dy,
  68.                 dx,
  69.                 cdy,
  70.                 lefty,
  71.                 righty;
  72.  
  73.     float      *m1,        /* four maps for bilinear interpolation     */
  74.                *m2,
  75.                *m3,
  76.                *m4;
  77.  
  78.     float     **maps;        /* pointers to intensity maps         */
  79.     float       equal_const;    /* constant to produce histogram equal. */
  80.  
  81.     int         i,
  82.                 j,
  83.                 k,
  84.                 loop;        /* utility infielders         */
  85.     int         x,
  86.                 y;        /* ditto             */
  87.  
  88.     register int *h;        /* fast pointer to histogram     */
  89.     register float *c;        /* fast pointer to cum. histo.   */
  90.     register float *m;        /* fast pointer to current     */
  91.                     /* intensity mapping         */
  92. #ifdef godot
  93.     extern int  report ();    /* status reporting routine     */
  94.     signal (SIGUSR1, report);    /* attach report routine to sig. */
  95. #endif
  96.  
  97. /* Preliminary setup    */
  98.     status = 0;
  99.     min = minmax[0];
  100.     max = minmax[1];
  101.  
  102. #ifdef DEBUG
  103.     fprintf (stderr,"ahe: min = %d\n", min);    /***/
  104.     fprintf (stderr,"ahe: max = %d\n", max);    /***/
  105. #endif
  106.  
  107.     maxlevels = max - min + 1;    /* runs from 0 to max + 1 */
  108.  
  109. #ifdef DEBUG
  110.      /*###*/ fprintf (stderr, "maxlevels = %d\n", maxlevels);
  111. #endif
  112.  
  113. /* calculate constants for matrix indexing    */
  114.     numline = dimv[2];        /* pixels in a line     */
  115.     regline = nreg[2];        /* regions in a line     */
  116.  
  117. /* Find the total number of regions to be considered:    */
  118.     numregions = 1;
  119.     for (i = 0; i < DIMNO; i++)
  120.     numregions *= nreg[i];
  121.  
  122. #ifdef DEBUG
  123.      /*###*/ fprintf (stderr, "maxlevels = %d\n", maxlevels) ;
  124.      /*###*/ fprintf (stderr, "numline = %d\n", numline);
  125.      /*###*/ fprintf (stderr, "regline = %d\n", regline);
  126.      /*###*/ fprintf (stderr, "numregions = %d\n", numregions);
  127. #endif
  128.  
  129. /* Allocate the histogram map.                    */
  130.     histo = (int *) Malloc(sizeof (int) * maxlevels);
  131.  
  132. /* Allocate the cummulative histogram map.            */
  133.     cumhisto = (float *) Malloc (sizeof (float) * maxlevels);
  134.  
  135. /* Allocate the intensity mapping tables:            */
  136.     maps = (float **) Malloc (sizeof (float *) * numregions);
  137.  
  138.     for (i = 0; i < numregions; i++)
  139.     *(maps + i) = (float *) Malloc (sizeof (float) * maxlevels);
  140.  
  141. /* 
  142.  * Compute the mappings for each region.  We step through
  143.  * the regions, calulating the histogram and the histogram
  144.  * equalizing mapping for each region.  The same histogram
  145.  * array is used over for each region, but each region has
  146.  * its own map. 
  147.  */
  148.  
  149. /* 
  150.  * Calculate region sizes.  If the image is not evenly divisible,
  151.  * the excess pixels are distributed among the first dimv mod n
  152.  * regions
  153.  */
  154.  
  155. /* Allocate the arrays to hold the size of each region:     */
  156.     for (i = 0; i < DIMNO; i++)
  157.     {
  158.     del[i] = (int *) Malloc (sizeof (int *) * nreg[i]);
  159.     curdel = del[i];
  160.     regsize = dimv[i] / nreg[i];
  161.     excess = dimv[i] % nreg[i];
  162.     for (j = 0; j < nreg[i]; j++)
  163.         *(curdel + j) = (j < excess) ? regsize + 1 : regsize;
  164.     }
  165.  
  166. /*
  167.  *    Determine cliplimit from range of intensities in image and
  168.  *    the region size.  Image size/(number intensity levels * number
  169.  *    of regions) gives approximation of number of pixels at each
  170.  *    intensity in a region.  Clipfrac is used to scale the result.
  171.  */
  172. #ifdef DEBUG
  173.     fprintf (stderr, "ahe: clipfrac %f\n", clipfrac);
  174. #endif
  175.     if (clipfrac > 0.)
  176.     {
  177.     cliplimit = clipfrac * (dimv[0] * dimv[1] * dimv[2]) /
  178.         (maxlevels * numregions) +.5;
  179. #ifdef DEBUG
  180.     fprintf (stderr, "ahe: cliplimit %d\n", cliplimit);
  181. #endif
  182.     }
  183.     else
  184.     /* 0 is flag for no clipping */
  185.     cliplimit = 0;
  186.  
  187. #ifdef DEBUG
  188.     fprintf (stderr, "ahecalc: clipfrac %f cliplimit %d\n", clipfrac, cliplimit);
  189. #endif
  190.  
  191. /* Set up to calculate the histogram and mappings.  */
  192.  
  193.     status = 0;
  194.     start[1] = start[2] = 0;
  195.     for (j = 0; j < nreg[1]; j++)    /* loop thru y regions     */
  196.     {
  197.         for (k = 0; k < nreg[2]; k++)    /* loop thru x regions     */
  198.         {
  199. #ifdef REPORT
  200.         printf ("status: %d/192", status);
  201.         fflush (stdout);
  202. #endif
  203.         sprintf (argv[0], "status:%d/192", status++);
  204.         newdel[1] = *(del[1] + j);    /* current region in y     */
  205.         newdel[2] = *(del[2] + k);    /* current region in x     */
  206.         equal_const = (double) (maxlevels - 1) / (newdel[1] * newdel[2]);
  207.  
  208. #ifdef DEBUG1
  209.          /*###*/ fprintf (stderr, " equal_const %f\n", equal_const);
  210.          /*###*/ fprintf (stderr, " newdel[1] %d\n", newdel[1]);
  211.          /*###*/ fprintf (stderr, " newdel[2] %d\n", newdel[2]);
  212. #endif
  213.  
  214.         /* zero histogram array                 */
  215.         for (h = histo; h < histo + maxlevels; h++)
  216.             *h = 0;
  217.         h = histo;
  218.  
  219.         /* calculate histogram for this region         */
  220.             for (y = start[1]; y < start[1] + newdel[1]; y++)
  221.             for (x = start[2]; x < start[2] + newdel[2]; x++)
  222.                 ++h[*(iptr1 (y, x))];
  223.  
  224.         /* clip histogram if user requested a limit         */
  225.         if (cliplimit > 0)
  226.         {
  227.             aheclip (histo, maxlevels, cliplimit, newdel[2], newdel[1],
  228.                  &ptemp, &p);
  229.  
  230.             /* calculate cumulative histogram             */
  231.             limit = cliplimit - ptemp;
  232.             if (*histo == limit)
  233.             *cumhisto = (float) cliplimit;
  234.             else
  235.             *cumhisto = (float) *histo + p;
  236.  
  237.             for (h = histo + 1, c = cumhisto + 1; h < histo + maxlevels;
  238.              h++, c++)
  239.             {
  240.             if (*h == limit)
  241.             {
  242.                 *c = *(c - 1) + (float) cliplimit;
  243.             }
  244.             else
  245.             {
  246.                 *c = *(c - 1) + (float) *h + p;
  247.             }
  248.             }
  249.         }
  250.         else
  251.         {
  252.             /* no clipping requested */
  253.             *cumhisto = (float) *h;
  254.             for (h = histo + 1, c = cumhisto + 1; h < histo + maxlevels;
  255.              h++, c++)
  256.             {
  257.             *c = *(c - 1) + (float) *h;
  258.  
  259.             }
  260.         }
  261.  
  262. #ifdef DEBUG1
  263.          /*###*/ fprintf (stderr, "j %d, k %d\n", j, k);
  264. #endif
  265.  
  266.         /* calculate intensity mapping number         */
  267.         m = mapval (j, k);
  268.  
  269.         *m = 0;
  270.         for (c = cumhisto; c < cumhisto + maxlevels; ++c, ++m)
  271.             *m = equal_const * (*c);
  272.  
  273.         /* increment starting point                 */
  274.         start[2] += newdel[2];
  275.         }
  276.         start[2] = 0;
  277.         start[1] += newdel[1];
  278.     }
  279.  
  280. /* Apply mappings to the image */
  281.     start[0] = 0;
  282.     loop = 0;
  283.  
  284.         start[1] = 0;
  285.         ymin = 0;
  286.  
  287.         for (j = -1; j < nreg[1]; j++)    /* y regions */
  288.         {
  289.         if (j == -1)
  290.         {
  291.             jval = 0;
  292.             dy = 0;
  293.             recy = 0;
  294.             ymax = (*(del[1])) / 2.0;
  295.         }
  296.         else
  297.         if (j == nreg[1] - 1)
  298.         {
  299.             jval = j - 1;
  300.             dy = 1;
  301.             recy = 0;
  302.             ymax = dimv[1];
  303.         }
  304.         else
  305.         {
  306.             start[1] += *(del[1] + j);
  307.             jval = j;
  308.             dy = 0;
  309.             recy = 1.0 / (float) (*(del[1] + j));
  310.             ymax = start[1] + (*(del[1] + j + 1)) / 2.0;
  311.         }
  312.  
  313.         for (y = ymin; y < ymax; y++)    /* y pixels in rgion */
  314.         {
  315. #ifdef REPORT
  316.             if ((loop % 4) == 0)
  317.             {
  318.             printf ("status: %d/192", status);
  319.             fflush (stdout);
  320.             }
  321.             if ((loop++ % 4) == 0)
  322.             sprintf (argv[0], "status:%d/192", status++);
  323. #endif
  324.             cdy = 1.0 - dy;
  325.             start[2] = 0;
  326.             xmin = 0; 
  327.  
  328.             for (k = -1; k < nreg[2]; k++)    /* x regions */
  329.             {
  330.             if (k == -1)
  331.             {
  332.                 kval = 0;
  333.                 dx = 0;
  334.                 recx = 0;
  335.                 xmax = (*(del[2])) / 2.0;
  336.             }
  337.             else
  338.             if (k == nreg[2] - 1)
  339.             {
  340.                 dx = 1;
  341.                 kval = k - 1;
  342.                 recx = 0;
  343.                 xmax = dimv[2];
  344.             }
  345.             else
  346.             {
  347.                 start[2] += *(del[2] + k);
  348.                 kval = k;
  349.                 dx = 0;
  350.                 recx = 1.0 / (float) (*(del[2] + k));
  351.                 xmax = start[2] + (*(del[2] + k + 1)) / 2.0;
  352.             }
  353.  
  354.             /* maps for the interpolation */
  355.             m1 = mapval (jval, kval);
  356.             m2 = mapval (jval, kval + 1);
  357.             m3 = mapval (jval + 1, kval);
  358.             m4 = mapval (jval + 1, kval + 1);
  359.  
  360.             i1 = iptr1 (y, xmin);
  361.             i2 = iptr2 (y, xmin);
  362. #ifdef FOO
  363.             fprintf (stderr,"x, y= %d %d ",xmin,y) ;
  364. #endif
  365.  
  366.             for (x = xmin; x < xmax; x++)    /* x pixels in region */
  367.             {
  368.                 lefty = dy * (m3[*i1]) + cdy * (m1[*i1]);
  369.  
  370.                 righty = dy * (m4[*i1]) + cdy * (m2[*i1]);
  371.  
  372.                 /* CHANGED */
  373.                 tmpint = (int) (righty * dx
  374.                      + lefty * (1.0 - dx));
  375.  
  376.                             if (tmpint > MAXCHAR)
  377.                 {
  378. #ifdef DEBUG
  379.                 fprintf (stderr,"%lf %lf %lf %lf %lf %d\n",lefty,righty,
  380.                      dx,dy,cdy,tmpint) ;
  381. #endif
  382.                                 tmpint = MAXCHAR ;
  383.                 }
  384.                             *i2 = (unsigned char) tmpint ;
  385.                 i1++;
  386.                 i2++;
  387.                 dx += recx;
  388.             }    /* end for x */
  389.             xmin = xmax;
  390.             }        /* end for k */
  391.             dy += recy;
  392.         }        /* end for y */
  393.         ymin = ymax;
  394.         }            /* end for j */
  395. }                /* end */
  396.  
  397. #ifdef godot
  398. report ()
  399. {
  400.     signal (SIGUSR1, SIG_IGN);
  401.     fprintf (stdout, "%3d 192\n", status);
  402.     fflush (stdout);
  403.     signal (SIGUSR1, report);
  404. }
  405. #endif
  406.